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Abstract. We attempt to study the characteristics of the different stellar populations present in the Galactic 
central region. A Monte Carlo method is used to simultaneously fit 11 thin disc and triaxial outer bulge density 
parameters on (K 3 , J-K s ) star count data in almost 100 windows from the DENIS near infrared large scale survey 
at -8°<1<12° and |b|<4°. Various bulge density profiles and luminosity functions were tested using a population 
synthesis scheme. The best models, selected by a maximum likelihood test, give the following description: the 
outer bulge is boxy, prolate, and oriented 10.6° ±3° with respect to the Sun - center direction. It seems that the 
main bulge population is not older than 10 Gyr, but this preliminary result needs further work to be confirmed. 
A significant central hole is found in the middle of the thin disc. We discuss these results in regard to previous 
findings and the scenario of bulge formation. 
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1. Introduction 

Due to the high extinction as well as the superimposition 
of different stellar populations in the region, the struc- 
ture of the inner Milky Way remains uncertain. In the 
last 10-15 years, the use of infrared data, less sensitive to 
extinction, made it possible to make much progress in our 
knowledge of the bulge region. However, some questions 
such as the precise orientation and the length of the outer 
bulge, as well as the existence and the length of the central 
disc hole have received contradictory answers. 

There is now a consensus that the outer bulge is triax- 
ial: Blitz & Spergel (1991) and Kent et al. (1991), using 2.4 
/im maps, and Nakada et al. (1991), using IRAS stars, de- 
tected asymmetries in longitude which they explained by 
the triaxiality of the outer bulge, also called the bar, with 
the near end at positive longitudes ; Binney et al. (1991) 
studied the gas kinematics and came to the same conclu- 
sion. Then, studies using COBE/DIRBE maps (e.g. Dwek 
et al. 1995; Binney et al. 1997; Freudenreich 1998; Lepine 
& Leroy 2000; Bissantz & Gerhard 2002), star counts (e.g. 
Nikolaev & Weinberg 1997; Stanek et al. 1997; Lopez- 
Corredoira et al. 2000), kinematics (e.g. Deguchi et al. 
2002; Zhao 1996) and microlensing (e.g. Zhao & Mao 1996) 
made it possible to deduce a more detailed description of 
the triaxial bulge/bar, and in particular gave estimations 
of its orientation in the Galactic plane. But, if all stud- 
ies converge to the description by a prolate shaped outer 



bulge with the major axis almost lying in the Galactic 
plane and the near end in the first quadrant, the esti- 
mated values of the angle of the major axis from the Sun 
direction vary between about 10° and 30°. Some stud- 
ies even present a bar with an angle around 40°-45° (e.g. 
Nakai 1992; Weinberg 1992; Deguchi et al. 1998; Sevenster 
et al. 1999; Lopez-Corredoira et al 2001). The reason for 
this difference is that they deal with two different things: 
Lopez-Corredoira et al. (2001), then Picaud et al. (2003), 
showed the existence of a structure at l=20°-27°, which 
may be at the top end of a long in-plane bar with an 
angle of about 40° from the Sun-center direction. This 
structure is distinct from the triaxial outer bulge studied 
in the present work, but the confusion often appears in 
the literature. 

Alard (2001) showed evidence for a nuclear bar in the 
most central parts of the Milky Way which may cor- 
respond to the inner bulge. But according to Ibata & 
Gilmore (1995), the inner bulge and the outer bulge may 
be two distinct populations. Therefore, in order to avoid 
any confusion with another bar or with the inner bulge, 
we will give the name outer bulge to the triaxial prolate 
structure observed in the inner region (|1| <10°) excluding 
the very central parts (|1| <1° and |b| <1°). 

Another important stellar population present in the 
outer bulge region is the thin disc. The shape of the inner 
thin disc is still a controversial issue. In particular, the 
existence of a truncation or a hole at its center is still in 
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doubt, even if there are more and more indications in this 
direction, in particular in external galaxies. For instance, 
Freeman (1970) classified the discs in two types, the first 
described by a single exponential, and the second modeled 
by the subtraction of two exponentials, i.e. with a central 
hole. Ohta et al. (1990), studying 6 early-type spiral galax- 
ies, showed that all had Freeman type II discs. According 
to Bagget et al. (1996), the proportion of galaxies hav- 
ing an inner truncated disc is at least twice as numerous 
in barred spirals as in non-barred ones. In our Galaxy, 
Freudenreich (1998) fitted his bulge and disc model on 
the COBE/DIRBE map and found a disc hole radius of 3 
kpc. More recently, Lepine & Leroy (2000) showed that a 
disc model including the central hole was more convenient 
to fit the brigthness distribution and the rotation curve, 
and Lopez-Corredoira et al. (2001) used a disc model with 
a truncation inside the Galactic ring at 3.7 kpc. 

In this paper, we have compared simulated star counts 
from the Besangon model of the Galaxy with DENIS 
near infrared data in almost one hundred windows of 
low extinction distributed between -12° and +8° in lon- 
gitude and -4° and +4° in latitude. Star counts are a 
better means than surface brightnesses to determine the 
3-dimensional density parameters: integrated fluxes are 
dominated by the brightest and closest stars while star 
count studies take into account a wider range of intrinsic 
stellar luminosities and distances. 

The article is organized as follows. In Sect. 2, we will 
describe the data (DENIS batches), the selection of low 
extinction windows and the selections made in color and 
magnitude. In Sect. 3, after a brief description of the 
Besangon model of the Galaxy, we will present the disc 
and outer bulge density distributions and luminosity func- 
tions (hereafter LF) used to compute the simulations and 
fit the parameters. In Sect. 4, we will describe the fitting 
method, based on Monte Carlo drawings and a maximum 
likelihood test. Results are given in Sect. 5 and compared 
with previous studies in Sect. 6. We conclude in Sect. 7. 

2. The Data 

2.1. DENIS Batches 

The DENIS (Deep Near Infrared Survey of Southern Sky) 
survey (Epchtein et al. 1997, Fouque et al. 2000) covers 
almost all the Southern Sky (97 % at the end of the survey) 
in strips of 30°xl2'. Photometric bands used are Gunn-I 
(0.85 ^m), J (1.25 ^m) and K s (2.15 ^m). 

In addition to the survey strips, specific observations 
called "batches" (Simon et al, in preparation) were made 
in 1998 in outer bulge and plane regions in smaller fields 
(about 2 deg 2 ). Data reductions were made at PDAC 
(Paris Observatory Data Analysis Center). A PSF fitting 
optimised for the crowded fields was used for the source 
extraction. All standard stars observed in a given night 
were taken for the determination of the photometric zero 
point, resulting in accuracies ranging for 0.08 at K s =8 to 
0.15 at K s =13, and 0.08 at J=10 to 0.15 at J=15. 



2.2. Windows of low extinction 

94 windows of about 15'xl5' between the longitudes [- 
12°; 8°] and the latitudes [-4°; 4°] were selected in DENIS 
batches. These windows were selected from the Schultheis 
et al. (1999) extinction map as having either a low or 
homogeneous extinction which is easy to model. For most 
windows, the extinction distribution along the line of sight 
was modeled using no diffuse extinction but 2 clouds (lo- 
calized extinction), one at 1 kpc (Sagittarius-Carina arm) 
and the other at about 4 kpc (Scutum-Centaurus arm). 
In a few windows, diffuse extinction was added. In each 
cloud, the Ay was estimated by comparing the quantilcs 
of J-K s color between data and simulations (see Sect. El). 
The reddening in each photometric band was taken from 
the extinction law of Mathis (1990). 

Unfortunately, as one can see in Fig. ^ there are al- 
most no windows selected very close to the Galactic plane 
(|b| <1°), because of the large extinction there. 

2.3. Photometric bands and star selections 

In the present study, only magnitude K s and color J-K s 
were used to compare observed star counts with simulated 
ones, the I band being too sensitive to extinction. 

The blue side of the color-magnitude diagram (here- 
after CMD) is mostly populated by foreground dwarfs, 
especially at faint magnitudes. Cuts were made in K s and 
J-Ks to reduce this contamination: we kept only stars be- 
tween 7.5 and an upper limit varying from field to field 
(12.5 on average) in K s , and with J-K s > J-K s /; im , J- 
Ks/zim being in the mean equal to 0.8 but also varying 
from field to field due to extinction. To ensure complete- 
ness, only stars below «15 mag in J have been kept. 

3. The model 

The Besangon model of the Galaxy has been used to com- 
pute simulations and compare them with the data. It is de- 
scribed in Sect. 13. ll The holed thin disc and triaxial outer 
bulge density models which have been fitted are presented 
there, as well as the luminosity function used. 

3.1. The Besancon model of the Galaxy 

The Besangon model of stellar population synthesis aims 
at giving a global 3-dimensional description of the Milky 
Way structure and evolution, including stellar populations 
such as thin disc, outer bulge, thick disc and spheroid, as 
well as dark halo and interstellar matter. Details can be 
found in Robin & Creze (1986), Bienayme et al. (1987) 
and Robin et al. (2003,2004). 

The approach of the Besangon model is semi-empirical: 
both theoretical schemes (stellar evolution, galactic evo- 
lution, galactic dynamics) and empirical constrained are 
used. Boltzmann and Poisson equations make it possi- 
ble to self-consistently constrain the disc scale height via 
the Galactic potential. Simulated catalogues of stars are 
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Fig. 1. Windows of low extinction (small filled boxes) and batches (long boxes). Grey levels of windows represent the 
Ay given by Schultheis et al. (1999). 
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built and give observables (apparent magnitudes, colors, 
proper motion, radial velocities) directly comparable with 
the data. Interstellar extinction, photometric errors and 
Poisson noise are also added to make simulations as close 
as possible to observations. 

The Besangon model has been used to determine the 
density laws and luminosity functions of the stellar pop- 
ulations: thin disc (Robin, Creze & Mohan 1992, Ruphy 
et al. 1996, Haywood et al. 1997), thick disc (Robin et al. 
1996, Reyle & Robin 2001), spheroid (Robin et al. 2000), 
and outer bulge in the present work. 

3.2. Holed Thin Disc 

The thin disc may contribute significantly to star counts in 
the Galactic central region. Hence its characteristics have 
to be fitted at the same time as the bulge parameters. 
The two disc shape parameters that have an effect on star 
counts of the inner region are the scale lengths of the disc 
(Rd) and of its central hole (Rh). Other parameters and 
structures such as outer radius, warp and flare are taken 
into account in the Besangon model but are not relevant 
here. 



is probably very patchy. Fitting its parameters would be 
difficult and inefficient. 

The old thin disc density distribution model follows 
the Einasto (1979) law: the distribution of each old disc 
component is described by an axisymmetric ellipsoid with 
an axis ratio depending on the age; the density law of 
the ellipsoid is described by the subtraction of 2 modified 
exponentials: 



,>., = ,>.,„ - [rxi>(- , 0.25 + (-|-) 2 )-exp(- JO.25 + (-^-) 2 ) 



with a 



R 2 



where: 



R and Z are the cylindrical galactocentric coordinates, 
e is the axis ratio of the ellipsoid. Table ^ gives the 
recently revised (Robin et al. 2003) axis ratios of the 
6 age components of the old thin disc. 
Rd is the scale length of the disc and is around 2.3-2.5 
kpc (Ruphy et al. 1996). 
Rh is the scale length of the hole. 
The normalization pd is deduced from the local lumi- 
nosity function ( Jahreifi et al, private communication) , 
assuming that the Sun is located at R0=8.5 kpc and 
Z =15 pc. Local densities are given in Tabled 



3.2.1. Density distribution 

The thin disc is divided into 7 age components: the first 
(age < 0.15 Gyr) is called young disc, the other six form- 
ing the old disc (0.15 Gyr < age < 10 Gyr). The young 
thin disc is not studied here, because its density is very 
low in comparison with the bulge and the old disc, and 



3.2.2. Luminosity Function 

A standard evolution model is used to produce the disc 
population, based on a set of evolutionary tracks, a con- 
stant Star Formation Rate (hereafter SFR) and a two- 
slope Initial Mass Function (hereafter IMF) (j){M) = 



4 



S. Picaud and A. C. Robin: The shape of the Galactic bulge from DENIS 



Table 1. Axis ratios and local densities of the six age 
components of the old thin disc. 



age (Gyr) 


e 


Pd (*.pc ' 5 ) 


0.15-1 


0.0268 


0.03146 


1-2 


0.0375 


0.02538 


2-3 


0.0551 


0.01887 


3-5 


0.0696 


0.02625 


5-7 


0.0785 


0.02037 


7-10 


0.0791 


0.02284 



Fig. 2. Luminosity functions of the thin disc in the K s 
band. The young disc corresponds to the dashed line, while 
the old one is represented by the dotted line. The solid 
line represents the total thin disc LF. On the ordinate: 
the decimal logarithms of the numbers of stars per 1 mag 
bin of absolute magnitude. 
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Absolute magnitude Ks 



A- M~ a with cv=1.6forM<l-M and a=3.0 for M>1M . 
The preliminary tuning of the disc evolution parame- 
ters against relevant observational data was described in 
Haywood et al. (1997) and further changes are explained 
in Robin et al. (2003). 

Magnitudes and colors in various filters are computed 
using semi-empirical model atmospheres from Lejeune et 
al. (1997, 1998). 

Fig. [21 presents the luminosity functions in K s of the 
young disc, the old disc and the thin disc in totality. This 
figure confirms how dominant the contribution of the old 
component on the thin disc is, in particular at the mag- 
nitudes [-5;-2] which are the most frequent absolute mag- 
nitudes in the simulations of the present study. Within 
the range of apparent magnitude and color used in the 
present study, observed stars have absolute K s magnitudes 
brighter than -1. 

3.3. Triaxial bulge 

To simulate the star counts, the bulge density law and a lu- 
minosity function have to be assumed. The different bulge 
density profiles used are presented in Sect. 13. 3.11 5 differ- 



ent luminosity functions have been tested, as explained in 
Sect. Eg 

3.3.1. Density distribution 
Orientation 

There is consensus that the bulge is triaxial. Three 
angles define the orientation: 

— 4>: orientation angle from the sun-center direction of 
the projection on the Galactic plane of the bulge major 
axis. 

— (3: pitch angle of the bulge major axis from the Galactic 
plane. In all previous studies, (3 was found to be very 
close to 0°. 

— 7: roll angle around the bulge major axis 

However, the third angle 7 is ill defined because the 
minor axes have similar scale lengths. Hence we prefered 
to fix it at 7=0°, and only (j> and f3 have been fitted. 

Density profiles 

Dwek et al. (1995) and Freudenreich (1998) fitted var- 
ious density distributions to the near infrared surface 
brightness observed using the Diffuse Infrared Background 
Experiment (DIRBE) of the Cosmic Background Explorer 
(COBE). The best-fitting models obtained at 2.2/im were 
the G2 and E\ functions by Dwek et al. (1995) and the 
S function in Freudenreich (1998). Stanek et al. (1997) 
tested the same functions of Dwek et al. (1995) using star 
counts of red clump giants. Their best-fitting model, con- 
sistent with the observed star counts as well as the surface 
brightnesses used by Dwek et al. (1995), was the E2 func- 
tion. 

We choose to test the 3 different functions used in these 
models, as presented in Table an exponential function 
(called E) like Ei and E2, a gaussian one (called G) like 
G2, and the S function of Freudenreich (1998), which is a 
sech 2 function. 

All the best density profiles obtained by Dwek et al. 
(1995), Stanek et al. (1997) and Freudenreich (1998) are 
included in the 3 functions described above: E\ corre- 
sponds to the E function with C11 = C± = 1 (dia- 
mond shape), E% is obtained using the E function and 
C\\=C±=2 (ellipsoidal shape), and the boxy profile G2 
using the G function and Cm =4 and C±=2. The best val- 
ues obtained in Freudenreich (1998), using the S function, 
are: C||=3.501±0.0016 and C_l=1.574±0.0014. 

The density function is then multiplied by the cut-off 
function f c (distances are given in kpc, and R c is called 
the cut-off radius): 



p = pxf c (R XY ) 1 with Rxy = VX 2 + Y 2 
Rxy <Rc=> MRxy) = 1 
Rxy >Rc=> MRxy) = exp (-( % y - fi ° ) 2 ) 
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Fig. 3. Definition of the angles to pass from the Galactic frame (x,y,z) to the bulge frame (X,Y,Z). The transformation 
consists of 2 consecutive rotations: the first is a clockwise rotation of <fi around the galactic vertical axis z, and the 
second is a clockwise rotation of (3 around the new axis y'. 
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Table 2. The 3 bulge density profiles used 




E p E = po x exp(—R s 
G pG — po x cxp(— 0.5 • Rl 
S ps = Po x sech 2 (— R a 



with f£ 11 = \\±\ c ± + |^| c ^1 c ii /c ^ + |^| c " 



+ Li. ~L + 

I J/0 I J 1 z I 



Finally, the two angles <j> and /?, the three scale lengths 
xq, yo, Zo, the density at the center po, the cut-off radius 
R c and the two coefficients Cm and C± are considered as 
free density parameters in the fitting process. 

3.3.2. Luminosity functions 

Star counts are a function of both the density law and 
the luminosity functions. Five theoretical bulge luminosity 
functions have been tested, all of them based on a Salpeter 
IMF (a— 2.35) and assuming a single epoch of formation 
(starburst) as well as a mean solar metallicity (Z«0.02). 
Only stellar evolutionary tracks and bulge age vary from 
one LF to another: 

— Three of them (Fig. 0J have been deduced from theo- 
retical isochrones by the Padova team (Girardi et al. 
2002), and three bulge ages were tested: 7.9 Gyr, 10 
Gyr and 12.6 Gyr, which we shall name Pad7.8, PadlO 
and Padl2.6 respectively. These luminosity functions 
were computed from the evolutionary tracks of Girardi 
et al. (2000) for the mass range O.15M - 7.0M Q 
and Bressan et al. (1993) for M > 7.OM , deduced 
from a combination of the Girardi et al. (2000) (with 
Z=0.019 and Y=0.273), and the Bertelli et al. (1994) 
(Z=0.02, Y=0.280) non-a-enhanced isochrones. Stellar 
atmosphere models were taken from the ATLAS9 the- 
oretical spectra library (Castelli et al. 1997). 

— The others two (Fig. EJ have been taken from evolu- 
tionary bulge synthesis models of Bruzual & Chariot 
(see Bruzual et al. 1997). Models were constrained on 
several bulge globular clusters. Two bulge ages were 
tested: 10 Gyr and 12 Gyr, respectively called BC10 
and BC12. Evolutionary tracks were deduced from 
the Padova 1994 set (Alongi et al. 1993; Bressan et 
al. 1993; Fagotto et al. 1994a,b; Girardi et al. 1996), 



with Z=0.02 and Y=0.280. Stellar atmosphere mod- 
els were taken from the semi-empirical Lejeune et al. 
(1997,1998) spectral library. 

While Bruzual & Chariot luminosity functions are de- 
duced from the Padova 1994 set of evolutionary tracks, 
those of Girardi et al. (2002) use the Padova 2000 models, 
which are a new version of the Padova 1994 isochrones. 
The model atmospheres also differ: Girardi et al. (2002) 
use the theoretical spectra from Castelli et al. (1997), 
while Bruzual & Chariot prefer the semi-empirical library 
from Lejeune et al. (1997,1998). 

Fig. HO compares Bruzual & Chariot and Girardi et 
al. (2002) (Padova) models for an age of 10 Gyr. In Fig. 
01 |5] and only dwarfs, subgiants, red giants, horizontal 
branch stars and AGB stars are shown. Planetary nebulae 
and white dwarfs, included in Bruzual & Chariot models, 
are unobservable or negligible in star counts in the range 
of apparent magnitude used in the present work. However, 
they will be taken into account in the calculation of the 
total number of bulge stars in Sect. 15.51 

The LFs differ only slightly while CMDs are sensi- 
tive to the assumed age. The main age effect is the po- 
sition of the turn-off, which corresponds to a color shift 
of about 0.05 mag for a change in age of 2 Gyr. The only 
non-negligible difference between Padova and Bruzual & 
Chariot CMDs is that the Padova asymptotic giant branch 
coincides with the Bruzual red giant branch. 

4. Fitting method 

The method used to fit old thin disc and triaxial bulge pa- 
rameters to DENIS data is based on two important points: 
firstly, comparisons between data and models are made 
with star counts in K s magnitude and J-K s color; sec- 
ondly, parameters are deduced from a fitting method using 
Monte Carlo drawings and maximum likelihood tests. 
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Fig. 4. Color-magnitude diagrams (left) and luminosity functions (right) from Padova models. The three models, Pad7.9, PadlO 
and Padl2.6, are shown in the graphs using dark grey, grey and light gray respectively. The luminosity functions are given in 
number of stars per 0.5 mag bin of absolute magnitude K s . 




Absolute colour J-Ks Absolute magnitude Ks 



Fig. 5. Color-magnitude diagrams (left) and luminosity functions (right) from the Bruzual & Chariot models. The darker grey 
corresponds to BC10 and the lighter grey to BC12. The luminosity functions are given in number of stars per 0.5 mag bin of 
absolute magnitude K s . 




Absolute colour J-Ks Absolute magnitude Ks 



4.1. Cuts and choice of bin size 

For each field, selected stars are distributed in 8x2 equally 
populated bins: 8 bins of magnitude, and 2 bins of color. 
We noticed that the noise relative to the least popu- 
lated bins increases and contributes too much to the 
global likelihood. This bias increases with the difference 
between data and model. To limit this bias, we group 
some close windows (always from the same batch) so that 
we have at least 70 stars in each bin of magnitude-color. 
Approximately 10% of the windows are grouped together 
by 2 or 3, and we obtain 88 groups or single windows at 
the end, for 94 windows in total. 

4.2. Monte Carlo drawings 

There are 11 parameters to fit: 

— Bulge orientation: (f>, (3 

— Bulge scale lengths: xq (major axis), j/o and zo 

— Bulge normalization po and cut-off radius R c 

— C|| and Cjl 

— Disc scale length Rd and hole scale length Rh 



An iterative scanning method to explore this 11- 
dimensional space of parameters would be too time con- 
suming. An alternative to save computer time would be 
to distribute the parameters in groups of 3 or 4 and fit 
them group by group. This method is still too slow and 
does not take into account correlations between parame- 
ters, and therefore can miss some solutions as well as pre- 
venting convergence in some cases. Hence we prefer to use 
another fitting method based on Monte Carlo drawings. 
This method was developed to solve non-linear equations 
of star kinematics (Oblak, 1983) and has been adapted for 
the present study. 

In the Monte Carlo method, parameters are drawn in 
the 11-dimensional parameter space. 

— At the first iteration, p points of the 11-D space (or 
sets of parameters) are drawn using uniform drawings. 
Ranges for uniform drawings are given in Table 

— At each next iteration, the m points having the high- 
est values of likelihood among all the points drawn 
since the beginning are extracted. Then p new points 
are determined using semi-gaussian drawings around 
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Fig. 6. Color-magnitude diagrams (left) and luminosity functions (right) from Padova (PadlO, dark grey) and Bruzual & 
Chariot (BC10, light grey) for an age of 10 Gyr. The luminosity functions are given in number of stars per 0.5 mag bin of 
absolute magnitude K s . 




the median of the m best points and along the axes 
denned by the eigenvectors of their correlation matrix. 
Formulae of median and semi-dispersions related to the 
semi-gaussian drawings are given in Appendix 1X1 
It can happen that drawn values go beyond restrictive 
limits. In this case, a new value is determined by a 
uniform drawing between the inferior/superior limit 
and the minimal/maximal value respectively reached 
by the m best points on the concerned parameter. More 
details about the limits are given in Sect. I4.2T1 

— The fit ends when there is no further progress in the 
likelihood convergence. A maximum of 20 iterations is 
admitted, but this limit is almost never reached. 

Various values of m and p have been tested and the 
values giving the best compromise between quality and 
rapidity of convergence were m=40 and p=300. 

4.2.1. Constraints on parameters 

Table |3 gives the ranges for the first iteration uniform 
drawings. These ranges were chosen to avoid unrealistic 
configurations, based on previous studies which agree on 
the following description: a prolate shaped bulge (small 
yo and Zq with respect to Xq) with the major axis closer 
to the Galactic plane (small 0) and its top end in the 
first quadrant (0° < <fi < 90°). The following three points, 
concerning the limits, should be noted: 

— In the case of the angle f3, the limits apply during 
the first iteration only. Subsequent iterations can go 
beyond it. 

— For certain parameters, such as the central density po 
and the bulge scale lengths xq, ya and zg, only the 
inferior limit cannot be passed, and is equal to + . 

— In the other cases (the parameters C|| and Cj_, the 
angle </>, the cut-off radius R c and the disc scale lengths 
Rd and R/j), the limits are kept for all iterations. 



4.2.2. Maximum likelihood test 

Maximum likelihood methods are better than minimum 
X 2 ones, especially in the case of small counts per bin, be- 
cause minimum x 2 methods assume that the distribution 
of observed star counts is a gaussian, and this hypothesis 
is false in such a case. Hence, we preferred to use the max- 
imum likehood to determine the best sets of parameters 
in the fitting method, even if normalized residuals and \ 2 
were also used as well as likelihood to evaluate, after the 
fits, the agreement between models and observations. 

The use of the maximum likelihood method only as- 
sumes that the likelihood is smooth enough and well de- 
fined around its local maxima, which is eventually the 
case. 

X 2 and likelihood formulas take into account the speci- 
ficities of the study, especially the presence of noise in the 
simulations. They are presented in Appendix 151 

4.3. Weighting 

Computing simulations for each set of parameters drawn 
and for all windows would be very time consuming. 
Rather, we start with initial simulations, then weightings 
of these simulations are applied to adapt them to different 
density laws. Initial simulations were done using the G2 
function (C||=4, C ± =2, 7=0°) from Dwek et al (1995) as 
the bulge density distribution, densities with E, G and S 
distributions being calculated by changing the weightings 
adequately. 

This method may lead to a bias if the values of param- 
eters used for the initial simulations are badly chosen: if 
for instance the initial configuration implies that there are 
no bulge stars in a field, then there will be no bulge con- 
tribution for star counts in this field even for a drawn set 
of parameters which forecasts the presence of the bulge in 
the field. This bias is avoided by making simulations with 
appropriate values of parameters, i.e. no disc hole (i?/ t =0 
kpc), bulge scale-lengths all equal to the upper limit of the 
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Table 3. Ranges for uniform drawings for the first iteration. Lower limits are given in the first line while upper limits 
are given in the second, and restrictions related to the limits applied in following iterations in the third, oo means no 
restriction, and [] corresponds to strict limits. 
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major scale- length xq (3 kpc), large bulge cut-off radius 
(5 kpc). 

A Poisson noise is added to the Besancon simulations 
to make them resemble the data as closely as possible. 
Initial simulations were computed with a density about 5 
times greater to increase their signal to noise ratio. 

4.4. Convergence 

4.4.1. Dispersions of best parameters 

Due to the Monte Carlo drawings, the m best sets of pa- 
rameters of one fit are not exactly equal but their values of 
likelihood are very close, or even the same. Furthermore, 
two different fits do not give the same mean of best val- 
ues. That is why 20 independent fits were made, the final 
result being the means and dispersions around the 20 xm 
best values. 

4.4.2. Quality of convergence 

A convergence test was made using simulations to estimate 
the correlations and quality of convergence and to identify 
possible biases in the procedure. Three simulations were 
produced to be used as data and the full procedure were 
applied to them to attempt to retrieve the parameters 
used for these simulations. The three sets of parameters 
used to make these simulations have been chosen from 
randomly drawn sets to give three sufficiently different 
configurations for the thin disc and the outer bulge. 

The main conclusions given by these three tests of con- 
vergence and correlations are the following: 

— The values of the parameters are generally retrieved 
with good precision, around 6% of the range of the 
initial drawings and 6% of the value , except for some 
parameters mentioned in the following lines. 

— The accuracy is good: of the 33 values of parameters to 
be found, a third are retrieved at less than 0.5 sigmas, 
about 60% at less than 1 sigma, about 80% at less than 
1.5, 88% at <2 and 97% at <2.5. 

— Disc scale lengths Rd and Rh are strongly anticorre- 
lated, which will be taken into account in the discus- 
sion. The disc parameters do not show any significant 
correlation with the bulge ones. However, Rd and Rh 
best values are sometimes a little biased with respect 
to the real values (varying from to 2.4 a difference). 



— Cm and Cj_ are slightly anti-correlated to yo and zo 
respectively. Moreover, Cm and C± are not precisely 
determined, having a dispersion between 0.5 and 0.7, 
and a difference with the original values of up to 1.6 
sigma. 

— (3 as well as the latitudes of the fields studied here being 
small, the line of sight is almost parallel to the plane, 
which implies a slight anti-correlation between xq and 
4>. Furthermore, <f> presents a slight correlation with 
the density at the center pq. Moreover, the precision 
of (j) depends on its value: it is good when <\> = 14°, less 
good when is small (8.5°) and much worse when <j> 
is high (37°). 

— The cut-off radius R c is badly constrained, with dis- 
persions varying between 0.3 and 0.85 kpc. 

5. Results of fits 

5.1. First round 

The first round of fits presents a great degeneracy: the 
angle 4> is found to converge towards two distinct solutions. 
While many fits converge to a median value in the range 
5°, 12.5°, many others give a value very close to 0°, which 
is inevitably biased because of the strict limit at 0° for 
this parameter. Medians and dispersions of the 2 groups 
of solutions are tabled. The cut between the two partitions 
has been fixed at 0=3.5°. 

The degeneracy does not depend on the tested density 
profile, but is more pronounced with Padova luminosity 
functions (especially for age 12.6 Gyr) than with Bruzual 
& Chariot LFs. 

5.1.1. Identification of badly fitted windows 

Fig. presents the map of the mean \r (square root of 
X 2 per bin ; see appendix IB. 21 for more explanations) by 
window associated to the best parameters for the Padova 
luminosity function at 10 Gyr. The maps related to the 
different LF are similar, and the badly fitted windows are 
mostly the same, even if the Xr values are a little smaller 
for the Bruzual & Chariot luminosity functions. Mean Xr 
are calculated over all the fits without taking into account 
the degeneracy. 

Two groups of bad fields can be extracted: 

— The first group of 8 windows is located at posi- 
tive longitudes (3° <1<6°) and negative latitudes (- 
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2° <b<0.5°), which contributes around 25% of the 
global likelihood. All these fields present a significant 
excess compared to model counts. A probable explana- 
tion is the following: these windows are amongst the 
most extincted of all the studied fields. Firstly, this 
implies that a poor modeling of the extinction distri- 
bution, fitted using the Bruzual & Chariot BC10 lumi- 
nosity function, has a larger effect on the star counts 
than expected. Secondly, because of this high extinc- 
tion, the observed stars are intrinsically brighter than 
in other directions, and the poor agreement may be 
caused by a problem in the model for this range of ab- 
solute magnitudes, which has almost no effect on the 
star counts in other locations. 
— A second group of 10 windows, located around the 
galactic center, contributes around 25% of the global 
likelihood. All these fields present a large deficiency 
in the star count compared to the model. This must 
be due to the absence of a stellar population in the 
Besancon model of the Galaxy, the inner bulge, dis- 
tinct from the outer bulge that we are studying here, 
and confined to the first 1 or 2 degrees around the 
Galactic center. The presence of such a central pop- 
ulation is mentioned by Ibata & Gilmore (1995) and 
Frogel et al. (1999). The nuclear bar evoked by Alard 
(2001) could correspond to the same population. 

These two different groups of badly fitted windows 
represent less than 20% of the total number of fields 
but contribute more than half of the global likelihood. 
Furthermore, the degeneracy in <j) is mostly implied by 



them. As the aim of this study is to obtain a large scale 
description of the outer bulge and thin disc populations, 
we prefer to remove these fields and make new fits without 
them. 

5.2. Second round 

Tables 0] to 03 show the best parameters related to the 
second round fits for the 5 luminosity functions. For a 
given LF and a given density profile, the medians and 
dispersions around these medians of the 20 xm best values 
of parameters have been calculated. 

One can see that the likelihood has decreased strongly. 
The degeneracy in (j) has almost disappeared with the 
Bruzual & Chariot BC10 and BC12 and Padova Pad7.9 
functions, and with the Padova PadlO LF when the S den- 
sity profile is used. It is less important but still present for 
the Padova PadlO (using E and G luminosity functions) 
and Padl2.6 functions. As the fits converge away from 0°, 
the limit between the two groups of the degeneracy has 
been moved from 0=3. 5° to 0=4.5°. 

The second group of fits (with <4.5°), which are 
much less numerous in the second round, is now considered 
an artefact, and only the other group of fits will be taken 
into account. 

5.3. Density profiles and luminosity functions 
Some conclusions can be deduced from Tables 0] to |H1 
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Table 4. Results of the second round of fits using the Pad7.9 luminosity function. The table is composed of 3 subtables, one 
for each density profile (E G and S) . Each subtable is formed by two couples of lines giving the median (/i) and dispersion (a) 
of the values. Each pair of lines corresponds to a group of solutions (see text): fits with <j> >4.5° (upper) and fits with <j) <4.5° 
(lower) . The number of fits of each group is given at the beginning of the first line. L and \ r mean global log-reduced likelihood 
and square root of \ 2 per Dm , respectively (see Appendix [B] for more explanations). 
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Table 5. Results of the second round of fits using the PadlO luminosity function. 
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Table 6. Results of the second round of fits using the Padl2.6 luminosity function. 
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— Density profiles: the gaussian function (G) (the least 
peaked at the center) is the worst of the three bulge 
density profiles tested. E and S functions reach similar 
conclusions but give slightly different best values, es- 
pecially for <j) and its correlated parameters. The best 
agreement is obtained with the S profile. 



— Luminosity functions: At similar ages, the Bruzual 
& Chariot luminosity functions give better agreement 
than the Padova (Girardi et al. 2002) ones. However, 
the best LF used come from Padova models with a 
bulge age of 7.9 Gyr. The Bruzual & Chariot isochrone 
at a similar age was not available. This is the youngest 
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Table 8. Results of the second round of fits usinj 


; the BC12 luminosity function. 
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tested value, which may mean that the bulge age might 
be smaller. 

— Influence of the LF and density profile on best param- 
eters: the thin disc scale length Rd and Rh determina- 
tions are robust over the different luminosity functions 
and density profiles which have been tested. This is 
not the same for the bulge parameters. For instance, 
the orientation angle <fi is different with the E profile 
and with the two others, and shows a significant de- 
pendency on the luminosity function: it is larger when 
the Bruzual & Chariot LFs are used, and decreases 
when the age increases. Correlated parameters such as 
R c or xo are also affected. 

5.4. Quality of fits and correlations 

Fig. |H1 shows the map of Xr (square root of the \ 2 per 
bin) associated to the fits using the Padova luminosity 
functions at 10 Gyr. The other maps related to other LFs 
are similar. One can see that the agreement is good (less 
than 3 sigmas) for all but one or two windows. 

5.4.1. Accuracy of the results 

The dispersions obtained for the majority of parameters 
are small, however some of them are not well constrained, 
such as the bulge shape coefficients C± and C\\ and the 
cut-off radius R c . Nevertheless, this is not a serious prob- 



lem. A small variation of these parameters does not change 
significantly the bulge spatial distribution: C\\ and C± give 
only a general indication of the shape of the outer bulge, 
and do not have significant influence on other parameters, 
except a little on the bulge minor scale lengths yo and zq 
with which they are correlated. As for the cut-off radius, 
the scale length xq being short, the density varies by 10% 
to 20% at the distance R c , and the effect of the cut-off is 
not very strong. 

5.4.2. Correlations 

Table El shows the matrix of correlations around the best 
parameters for the second round of fits (group with <f> > 
4.5° only). The choice of the used luminosity function and 
density profile does not have any significant influence on 
the correlations, the matrix corresponds to the mean val- 
ues over all the (LF, profile) pairs. For a given LF and a 
given density profile, correlations are computed, using all 
the fits, as follows: let (£ fe ) be the mx 20 best points and Wk 
their weight (deduced from their likelihood) ; let us take 
two axis i and j, the coordinates & and £j of the points 
on these two axes, their weighted means m, = E/c Wk ■ 
and ruj — J2k w k ' £f> aim their weighted dispersions 

Si = Efc^fc ' (£i - m i) 2 and s 3 = Efe^fc ' (£jf " m 3? 
; then the correlation between the parameters i and j is 

equal to — . 
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Fig. 8. Map of the mean \ r by window related to the PadlO luminosity function. 
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Table 9. Mean values over the 5 luminosity functions and the 3 density profiles of the correlations of the first group fits of the 
second round. Values over ±0.7 are written in boldface. 
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The main conclusions given by these correlations are 
the following: 

— The disc scale lengths Rd and Rh are not correlated 
with bulge parameters but are strongly anticorrelated: 
the variation in the density due to decreasing one of 
the two parameters can be compensated by increas- 
ing the other one. This results in an incertainty of the 
hole scale length Rh, though without questioning the 
existence of the hole. 

— As we have already shown with the tests (see Sect. 
14.4.2(1 . the orientation angle <fi is anti-correlated with 
parameters such as xq and yo. 

— The bulge major scale length xq and central density 
P0) whose variations have opposite effects on the num- 
ber of stars at a given distance, are also strongly anti- 
correlated. 



5.5. Best parameters 

Table ITU1 gives values of the fitted parameters for the best 
model (using the Freudenreich (1998) sech 2 density profile 
and the Girardi et al. 2002 (Padova) luminosity function 
with an age of 7.9 Gyr), as well as the ones obtained by 
calculating the mean of the best values over all the pairs of 
density profile / luminosity function. These two sets of pa- 
rameters are consistent with each other. We can therefore 
deduce that the configuration described here is robust, and 
does not depend on the choice of the LF and the density 
profile. 

This configuration is the following: 

— Outer bulge: the outer bulge is prolate, with axis ratios 
1 : 0.31±0.4 : 0.26±0.3 (values from mean best mod- 
els), and seems to be boxy in all directions (Cj_ >2, 
C|| >2). Its major axis lies almost in the Galactic plane 
(/? is consistent with the value 0°) and is oriented about 
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Table 10. Best values obtained from the second round of fits. The first double line corresponds to the model giving the best 
agreement with data, i.e. with the S density profile and the Pad7.9 luminosity function. The values of the second one are the 
mean and dispersions of the best set of parameters obtained with each pair density profile / luminosity function. Meaning of the 
parameters: <f> gives the orientation of the bulge major axis with respect to the Sun - center direction ; /3 is the angle between 
this bulge major axis and the Galactic plane ; the bulge major scale length xo, whose significance depends on the density profile, 
has been replaced by xo, which is the distance on the major axis at which the density is equal to 38.6% of the central one ; 
r y = and r z = |^ are respectively the axis ratios associated to the bulge minor scale lengths yo and zq ; iVtot gives the total 
number of bulge stars ; Rc is the cut-off radius of the outer bulge, and Cy and C± correspond to its shape coefficients ; R d and 
Rh are respectively the scale lengths of the thin disc and of its central hole. 
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10° with respect to the Sun - center direction. It con- 
tains 82-10 9 ±20-10 9 stars, which corresponds to a mass 
of 2.4±O.6-1O 1O M . 
— Thin disc: the thin disc shows a large hole at its center. 
With a disc scale length of about 2.4 kpc and a hole 
scale length of about 1.3 kpc, the in-plane density goes 
from zero at the center to its maximum at a distance 
of 2 kpc, and then decreases until the cut-off at 14 kpc. 
However, the behaviour of the disc in the inner 100 pc 
is not constrained here. 

6. Discussion 

Many studies have been made of the structure of the outer 
bulge region. Here, we compare the description obtained 
by our fits with those found in the literature. 

6.1. Disc hole 

The strong anti-correlation between the two disc param- 
eters Rd and Rh forces us to be careful with the results, 
which may be biased. Nevertheless, the hole radius value 
(assumed to be the radius of the maximal disc density), ~2 
kpc, is large enough to be used as evidence of the presence 
of a central hole in the inner thin disc. 

This conclusion is in agreement with the studies of 
external galaxies by Ohta et al. (1990) and Baggct ct al. 
(1996), who argued for the existence of holed or inner trun- 
cated discs in most spiral barred galaxies, as well as the 
Galactic disc models of Freudenreich (1998) who studied 
the CORBE/DIRBE map and obtained a hole radius of 3 
kpc, Lopez-Corredoira et al. (2001) who analysed DENIS 
data and found a maximum disc density at about 2 kpc 
from the Galactic center, and Lepine & Leroy (2000) who 
found that a model including a central hole is in better 
agreement with the COBE/DIRBE surface brightness dis- 
tribution and the rotation curve. 

However, a central hole or an inner truncation are not 
the only things which can give a smaller disc density in 
the plane close to the center. An alternative explanation 
is a vertical flare in the inner disc. The scale height of the 



disc may be higher there as an effect of the bar potential. 
Therefore, assuming that the surface density of the disc 
is constant, this flare leads to a decreasing of the disc 
density in and close to the Galactic plane, as does a central 
hole. For instance, Lopez-Corredoira et al. (2004) found 
similar agreements with their data in the plane in the 
range 2-8 kpc using either a flaring disc or a holed one. As 
a flaring disc model does not change the global density, 
contrary to a holed disc model, a further study of the 
vertical distribution of the inner thin disc might make it 
possible to settle between these two alternative models. 

6.2. Outer bulge age 

All the tested bulge luminosity functions have been built 
using the same scheme: the bulge is composed of only one 
generation of stars, which are rather old, and the mean 
metallicity has been assumed solar but with a large dis- 
persion of 0.5 dex. These hypotheses are consistent with 
most constraints found in the literature. 

The best agreement in our fits is obtained with the 
luminosity function from Girardi et al. (2002) with an age 
of 7.9 Gyr. This is the youngest tested age. Therefore, we 
can assume that the bulge age is at least younger than 10 
Gyr, and perhaps even younger than 8 Gyr, which is in 
contradiction with the results of Zoccali et al. (2003) who 
propose 10 Gyr as a minimum, but is in favour of Cole 
& Weinberg (2002) who claimed that the outer bulge is 
not older than 6 Gyr. New luminosity functions with a 
younger age from Bruzual & Chariot (2003) as well as 
from Girardi et al. (2002) (Padova) are now available and 
we shall attempt to use them in the future to test the 
hypothesis of a younger age, or several mixed generations 
of stars, as well as the influence of the assumed bulge 
metallicity. 

6.3. Outer bulge shape 

Our best density profiles are the sech 2 function from 
Freudenreich (1998) and the exponential one proposed by 
Stanck ct al. (1997), with a preference for the sech 2 pro- 
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file. On the contrary, the gaussian function, the best model 
found by Dwek et al. (1995), gives the worst agreements. 

The outer bulge as described by our best parameters 
is very boxy. This is consistent with most other works on 
the subject, for instance Weiland et al. (1994) and later 
articles (Dwek et al. 1995, Freudenreich 1998 ...) from 
the same authors, which studied integrated luminosities 
from the COBE/DIRBE near infrared map. Moreover, 
the axis ratios, 1:0.30±0.02:0.25±0.01, which show that 
the triaxial bulge is very prolate, are similar to values ob- 
tained by Freudenreich (1998) (1:0.37:0.26), Dwek et al. 
(1995), Bissantz & Gerhard (2002) (1:0.3-0.4:0.3, also from 
the COBE/DIRBE map) and Weiner & Sellwood (1999) 
(1:0.33:0.33, using HI and CO data). This description of 
the outer bulge as a boxy prolate spheroid makes it similar 
to a bar. 

Concerning its half-length, the cut-off radius R c is fit- 
ted to 3.71±0.71 kpc. Using #0=1.82 kpc, assuming no 
cut-off, the bulge density on the major axis is equal to 3% 
of the central density at 4.4 kpc (R c +la) from the center, 
and 14% at 3 kpc (R c -la). This means that a cut-off at 
3.71 kpc is too far to be very pronounced. This explains 
the small precision in R c . The existence of such a cut-off 
may be related to the disc dynamics. A cutoff may appear 
at the corotation associated with the spiral arms or the 
molecular ring. However, the position found here is a bit 
too far out and inaccurate to allow us to conclude that it 
is a dynamical cut-off. 

6.4. Outer bulge orientation 

As in most of the other related works, the angle found 
with our fits is very small and consistent with 0°, which 
means that the outer bulge major axis almost lies in the 
Galactic plane. 

On the contrary, our estimation of <f> (i.e. the angle be- 
tween the bulge major axis and the Sun - center direction), 
10.6°±3°, differs from some previous works. These stud- 
ies, based on COBE/DIRBE surface brightness (Dwek et 
al. 1995, Binney et al. 1997, Bissantz & Gerhard 2002), 
kinematics studies (Feast & Whitelock 2000, Bissantz et 
al. 2002), or IRAS (Nikolaev & Weinberg 1997, Deguchi 
et al. 2002) and OGLE star counts analyses (Stanek et al. 
1997), give values around 20°. However our estimation is 
consistent with the 12°±6° obtained by Lopez-Corredoira 
et al. (2000), who analysed star counts from near infrared 
large scale survey as we did, and compatible with the 14° 
found by Freudenreich (1998) and Lepine & Leroy(200), 
both using the COBE/DIRBE map, or even the 16°±2° 
obtained by Binney et al. (1991) from gas kinematics. 

It is not clear why different studies coming from the 
same data (COBE/DIRBE) give somewhat different re- 
sults, which vary between (f> = 20°-30° and 4> = 14°-16°. 
The low sensitivity to bulge stars at negative longitudes, 
because of their large distance from us, may explain the 
difficulty in determining precisely the bulge orientation 
angle. Furthermore, it should be noted that these studies 



were based on integrated flux density, hence they were less 
sensitive to the distribution of stars along the line of sight 
than the present star count analysis. 

7. Conclusion 

We constructed a Monte Carlo fitting method to deter- 
mine the spatial distribution of outer bulge and disc stars 
from comparisons between DENIS K s and J-K s stars 
counts and simulations from the Besancon model of the 
Galaxy. Our best parameters show a thin disc that has a 
central hole, as often seen in barred spiral galaxies, and a 
boxy prolate outer bulge with a major axis lying almost in 
the Galactic plane, as obtained by most of the other stud- 
ies. However, our orientation angle (with respect to the 
Sun - center direction), <j> — 10.6° ± 3°, is slightly smaller 
than the mean value found in the literature. 

The best fit among the different bulge ages which have 
been tested is the youngest: 7.9 Gyr. Future fits involving 
deeper data and more fields close to the Galactic plane are 
planned to test other luminosity functions, with younger 
ages and varying metallicities. 

This study of the inner Galaxy is based on near in- 
frared data, using low extinction windows. This approach 
does not allow us to have fields very close to the Galactic 
plane and center. The use of mid-infrared observations, 
like ISOGAL (Omont et al. 2003) or GLIMPSE (Benjamin 
et al. 2003), or deep near infrared surveys like WIRCAM 
at the CFHT, soon to be available, will allow us to reach 
these fields. It would be interesting to be able to include 
these fields, especially to study the population of the inner 
bulge. 

The boxy prolate shape and a possible young bulge age 
are compatible with a bar structure, and may be explained 
by a reduced time scale for stellar formation. However this 
photometric study is not sufficient to settle the nature 
of the formation scenario of the outer bulge. Kinematical 
data would be needed to confirm the bar nature of the 
outer bulge, and thence its scenario of formation. 

Appendix A: Semi-gaussian drawings 

At each iteration of the fitting program, gaussian drawings 
determining the new p points were made using different 
values of dispersion at the right and at the left of the 
median point on the eigenframe axes. The aim was to de- 
scribe as fully as possible the morphology of the maximum 
likelihood region and reproduce it in the new drawings. 

Let (£ fe , k=l..m) be a family of m points of the nor- 
malized 11-D space of parameters. These points are sorted 
with respect to their reduced log-likclihood (see Appendix 
IbTi . L\ and L m being respectively the maximal and mini- 
mal likelihoods, and weighted, the weight Wk being defined 
by: 

w k = exp(-1.5- —J 

Li m — Li l 

This formula allows us to enhance (under a certain limit) 
the contribution of best points in the following formulae. 
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Let £ be the median of the m weighted points and Vj be 
the eigenvector j of their covariance matrix. Their values 
on the i coordinate in the initial frame are respectively 
and Vij . We call £^ the coordinate j of £ k in the eigenframe 

(?,{T^J=i..n}). 

We define the semi-dispersion o~ and <r+ by: 



Let (Q ) Z=l..jj,j=l..ll) be a family of px 11 coordinates 
obtained using a gaussian drawing of mean and disper- 
sion 1. The coordinates (£j,Z=l..p) in the initial frame of 
the p new points to be tested are then deduced from the 
following formula: 

£ = ^ + ct -E^+t+E^ 

M<0 % j>0 



One has: f{Zi/z t ) = ^f{Z Qi /z Ql ). According to Bayes' 
Theorem: 



f(Z i/z 0i ) = 



P(zo l /Zo i )f(Z 0i ) 
f P(z 0i /Z 0i )f(Z 0i )dZ 0i 



(B.3) 



Having no information a priori on Zq^ one can say 
that all f(Z 0i ) are equal 2 (at least around Zot, where 
P(zoi/ Zoi) is not negligible), so one can simplify Eq. 
El 

P{zoi/Z 0i ) 



f{Zoi/ Zoi) 



P(z 0l /Z 0i 
so: 

f(Zi/zi) -- 



J P( z 0i/ Zoi)dZ 0i 
LjZo^-e^, J P(z 0i /Z 0i )dZ 0i = 1, 



1 1 

a z 0i 



;Zoi 



(B.4) 



On substituting for P{y l /Z i ) from Eq.ESJ for f(Zi/zi) 
from Ea. IB. 41 putting these into Ea. IB. II and replacing 
Zi by aZoi, one finds (with £ = (1 + a)Z oi ): 



AaZ 0l )y>e- 



+i 



y z lz 0l \ (l + a)w+*o, 
{Vi + zoj)\ a Vi 

y t \z Qz l (1 + a)*+*>n+ 1 
{Vi + ZQi)\ ZOi+1 z\ 

~ z 0. 



a z 0i < 



Appendix B: Likelihood and \ 2 formulae 

Generally, when the agreement between observations and A 
an analytic model is tested, only the data has Poisson 
noise, and usual likelihood or \ 2 formulae are used. In 
the case of the Besangon model, initial simulations also 
have a Poisson noise. This particularity must be taken 
into account in the formulas of likelihood and \ 2 j as we ll 
as the fact that the model counts are deduced from initial 
simulations using weightings. 

Let i identify one of the bins. Correlations between ad- 
jacent bins are neglected. Let yi represent data star counts 

and Zi be model counts, yi (integer) follows a Poisson law Actually, we usually use the reduced log-likclihood, 
around Yi (unknown) . Zi (real) is obtained from an initial which is the natural logarithm of the likelihood minus the 
simulation star count Zi = azQ i: a being the weight 1 and likelihood computed with Zi replaced by yc 
ZQi (integer) following a Poisson law around Zoi = —Z^ 
{Zi and Z 0i being real). 



vJ-zqJ. 



( Zi + z 0i )v^h+i 



In 



{Vi 



B.l. Reduced log-likelihood 

Hereafter, probabilities (for discrete variables) will be 
noted P and densities of probabilities (for continuous vari- 
ables) will be noted /. 

By definition, the likelihood Ci (Kendall & Stuart 
1973) is the probability for an observed star count to 
take the value yi, assuming that the model is correct, i.e. 
Yi=Zi=aZoi- One has: 



Vi [ -Zoi { - 



Z 0iV- z 0i 
z 0i 



, {yt + z 0i y. ZOi 

m ; —z 0i u ' 



Vi 



Vi [ -Zoi { - 



(yi + z Ql )y-+^+ 1 

zoi + Zi 



= Vi ln( — ) - (Vi + zoi + 1) ln( 

Vi z 



Vi 



Ci = P{yi/zi) 



P{y i IZ i )f{Z i jz i )dZ i 



(B.l) 



Finally, by summing all the bins, one obtains: 



L = £[ ? a In(-) - (Vi + *K + 1) ln(^i±*)] 
^ y% z 0i + yi 



yi follows a Poisson law around Zi, so: 



1 



P(Vi/Zi) = —Z^e- Z * 



(B.2) 



The relative difference between the reduced log- 
likelihood obtained here and the usual one (i.e. without 
model noise: Lq — YliiUi m (™) — ( Zi ~~ Vi)) var i es between 
8% and 12%, which is not negligible. 



1 The case a — is excluded. 



2 Bayes' Postulate 
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B.2. X r 

yi is a Poisson random variable with an average and vari- 
ance Yi. This variance is unknown and its only approxi- 
mation available is the value yi. 

Zi is a random variable with an average aZ^i and a 
z 2 

variance a Z Qi = This variance is unknown and its 

z 2 

only approximation available is the value 

? z 2 

So, yi — Zi is a variable with a variance 07 s» ^ + 
One therefore obtains the following formula of normalized 
residuals: 

Vi - Zi 
r, = — 




\ " the end, by summing over all the color magnitude 
bins of a window, one obtains the x 2 per field (or group 
of fields): 



Fig. B.l. Reduced log-likelihoods (ordinates) versus 
square root of x 2 per bin. 



o 



o 
o 



o 



o 

X 

Lfj 




12 3 



= E^ 2 = E 



{Vi - Zjf 
1 z i 



However, we would rather use the square root of the 



X 2 per bin Xrf — y ^- (with n = number of bins) wich 
corresponds to the mean dispersion needed for random 
variables centered on the yi and used as model counts to 
obtain on average the same value of x 2 ■ 

By summing over all the fields or groups of fields 
(n/ =88 in total), one obtains the global x 2 and the square 
root of the x : 

Xr = ^.Xrf 

B.3. Reduced log-likelihoods vs \r 

If the likelihood has been preferred to Xr to extract best 
sets of parameters in the fitting method, the Xr is a more 
practical tool than likelihood to 'read' the quality of a fit, 
because it gives directly, in number of sigma, the mean 
distance between model and data. That is why the values 
of Xr have been calculated as well as likelihoods. 

Fig. lB.ll oresents the values of reduced log-likelihoods 
L versus the square root of x 2 P er bin Xr obtained by 
random drawings. 
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